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ABSTRACT 

A simphfied model of natural convection, similar to the Lorenz (1963) system, is com- 
pared to computational fluid dynamics simulations of a thermosyphon in order to test 
data assimilation methods and better understand the dynamics of convection. The 
thermosyphon is represented by a long time flow simulation, which serves as a refer- 
ence "truth" . Forecasts are then made using the Lorenz-like model and synchronized 
to noisy and limited observations of the truth using data assimilation. The result- 
ing analysis is observed to infer dynamics absent from the model when using short 
assimilation windows. 

Furthermore, chaotic flow reversal occurrence and residency times in each ro- 
tational state are forecast using analysis data. Flow reversals have been successfully 
forecast in the related Lorenz system, as part of a perfect model experiment, but never 
in the presence of significant model error or unobserved variables. Finally, we provide 
new details concerning the fluid dynamical processes present in the thermosyphon 
during these flow reversals. 



1 Introduction 



Forecasting methodologies, traditionally motivated by 
numerical weather prediction (NWP), can find applica- 
tions in other fields such as engi n eering ( Savelv et al.l . 



1972fl . finance (ISornette and Zhoul 120061 : iBollen et all 



2011 ), epidemiology ( Ginsberg et al.l . I2OO9I ). and marketing 
( Asur and Huberman . I2OIOI ). Techniques borrowed from the 
weather forecasting community may prove to be powerful for 
forecasting these other types of complex systems. Fluid sys- 
tems can be particularly challenging due to dynamics taking 
place at multiple interacting spatial and temporal scales. 
However, because of their relationship to NWP, fluid sys- 
tems are among the most studied in the context of forecast- 
ing. 

In this paper, we show that the flow in a computa- 
tional fluid dynamics (CFD) simulated thermosyphon un- 
dergoing chaotic convection can be accurately forecast using 
an ord inary different ial equation (ODE) model akin to the 
classic iLorenj (|l963l ) system. The thermosyphon, a type of 
natural convection loop or non-mechanical heat pump, can 
be likened to a toy model of climate. Thermosyphons are 
used in solar water heaters (|Belessiotis and Mathioulakisl . 
20021 ). cooling systems for computers l|Beitelmal and Patell . 
2002!) . roads and railways that c ross permafrost (iLustgartenl . 



2006 



nuclear power plants (iDetman and Whipd . 



Beine et al.l. Il992l : iKwant and Boardmanl . Il992l) . and other 
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industrial applications. In these heat pumps, buoyant forces 
move fluid through a closed loop, and at high amounts of 
forcing they can exhibit com plex aperiodic behavior. As first 
sugg ested bv lLorenj (|l963l ). this is illustrative of the unpre- 
dictable convection behavior observed in weather and cli- 
mate dynamics. 

Synthetic observations of the thermosyphon are com- 
bined with model data to produce new forecasts in the 
process known as data assimilation (DA). DA is a generic 
method of combining observations with past forecasts to 
produce the analysis, an approximately optimal initial con- 
dition (IC) for the next forecast cycle. Another interpreta- 
tion of the analysis is that it is a "best guess" for the true 
system state as represented in the phase space of the model. 
DA can be used as a platform for the reanalysis of past ob- 
servations, in which the dynamical model plays a key role 
in constraining the s tate estimates to be physically realistic 
(|Compo et al.i. l2006h . 

In the present study, we use an analysis of simulated 
thermosyphon mass fiow rate data to explain the heat trans- 
port processes occurring during chaotic flow reversals, and 
to inform empirical forecasts of the occurrence of these flow 
reversals. Although the flow reversals are chaotic, we show 
they have short-term predictability, quantifying the extent 
to which this is possible with our methods. 

This paper is structured in the following way: In Sec. [51 
we explain the CFD simulation used to generate a synthetic 
true state or "nature run" of the thermosyphon and the sep- 
arate forecasting model. In Sec. [3] we present an overview of 
how DA was applied to this experiment and its performance. 
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Figure 1. The thermosyphon has a simple circular geometry. The 
bottom wall is heated to a constant hot temperature Tf^ while the 
top wall is maintained at the temperature Tc, creating a temper- 
ature inversion of hot fluid below cold fluid. If conduction alone 
cannot stabilize this temperature inversion, then the fluid will 
begin to rotate and convection becomes the dominant process of 
heat transfer. The relevant model state variables are proportional 
to the bulk fluid velocity u and the temperature difference across 
the loop ATa—g. For counterclockwise flow, as indicated by the 
arrow near 9 o'clock, fluid velocity u > and AT3_9 is typically 
> 0. The radius ratio R/r = 24 used in our experiments is shown 
to scale. 



In Sec.|4]we explain and present the results for flow reversal 
and rotational state residency time forecasts. Finally, Sec. [5] 
contains concluding remarks. In Appendices Sl-4 in the Sup- 
porting Information we present a derivation of the model, 
detail the tuning of model parameters, and explain in detail 
the DA methods used. 



2 Models and the data assimilation algorithm 

Following previo us ex periments that examined 
the pe riodic (iKelleij 1 19661) a nd chaotic (IWelandeij . 
1967 : ICreveling et al. . 19751: [Gorman and Widmann, 
1984 : iGorman et all Il986l : Ehrhard and Mulior 



20061 : 



Ridouane et al 



1990l: lYuen and Bad . Il999l: Ijiang and Shoiil. 



2003 



Burroughs et al. . 200E ; iDesravaud et al.l . l2006l : lYang et al. 



20091 ) 



behavior of toroidal ther- 
we also consider a circular thermosyphon 



mosyphons 

geometry. Picture a vertically-oriented hula hoop, as shown 
in Fig. [1] An imposed wall temperature Th on the lower half 
of the loop (— -I < (j> < ^) heats the fluid contained in this 
section. Similarly, a wall temperature Tc < Th is imposed 
on the upper half (|- < (j) < ^)to cool the upper section 
(Fig. [1]). The forcing, proportional to the temperature 
difference AT™ = Th — Tc, is constant. We focus on the case 
of developed flow, ignoring transient behavior. 

The behavior of the fluid can be qualitatively under- 
stood as follows. As the heating parameter is increased, the 
flow behavior transitions from a conduction state (conduct- 



ing equilibrium) to a steady, unidirectional state of convec- 
tion (convecting equilibrium). No particular rotational state 
(clockwise, CW, or counterclockwise, CCW) is favored due 
to symmetry. At still higher heating values, chaotic flow os- 
cillations can be observed. In the chaotic regime, the flow is 
observed to oscillate around one unstable convecting equi- 
librium state until flow reversal. Each flow reversals causes 
the system to transition between CW and CCW rotational 
states. 



2.1 Thermosyphon simulation 

The reference state of the thermosyphon is represented 
by a CFD-based numerical simulation in two spatial di- 
mensions (2D). The details of the computational model 
have been des c ribed in detail in a previous study by 
iRidouane et"aLl (|2009l ): however, for completeness, we sum- 
marize here its essential elements. 

It is assumed that the temperature differential AT„ is 
sufficiently small so that temperature-dependent variations 
of material properties can be regarded as negligible, save for 
the density. The standard Boussinesq approximation is in- 
voked and all ffuid properties are assumed to be constant 
and evaluated at the reference temperature {Th + Tc)/2. 
The ffow is assumed to be laminar, two-dimensional, with 
negligible viscous dissipation due to low velocities. Under 
these circumstances, the governing dimensionless equations 
are the unsteady, 2D laminar Navier-Stokes equations along 
with the energy equation and equation of state for the den- 
sity. No slip velocity boundary conditions are imposed on 
the walls and isothermal boundary conditions of Th and Tc 
are imposed on the heated and cooled lower and upper walls, 
respectively. 

The dimensionless control parameter for convection is 
the Rayleigh number, defined here as 



Ra : 



Sg-yr^'ATu, 



(1) 



where g is the gravitational acceleration, 7 is the thermal 
expansion coefficient, v is the kinematic viscosity, and k is 
the thermal diffusivity. 

The one dimensionless geometric parameter is the ratio 
of major (loop) radius R to minor (tube) radius r, hereafter 
referred to as the radius ratio. Consistent with the previous 
study, the dimensions of the loop are chosen with i? = 36 
cm and r = 1.5 cm to yield a radius ratio of 24. 

As in the classic Rayleigh-Benard problem, the Rayleigh 
number determines the onset of convection in the ther- 
mosyphon. For the numerical simulations on this fixed ge- 
ometry, a range of Rayleigh numbers can be imposed by 
varying the value of the gravitational acceleration. As the 
Rayleigh number is increased from zero, the flow behavior 
transitions from a stationary, conduction state to a steady, 
unidirectional state of convection. At still higher values of 
Ra, chaotic fiow oscillations can be observed. Unless other- 
wise indicated, the simulation results presented in this paper 
correspond to a value of Ra= 1.5 x 10'', which is within the 
chaotic regime. 

All numerical simulat ions were performed using the 
commercial CFD software IaNSYS Fluent! (|2Q06l ). which is 
based on the finite- volume method. (An example of the out- 
put is shown in Fig. [8l in the discussion of fiow reversals.) 
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During the course of the simulations, the time- varying mass 
flow rate, a scalar denoted by q and proportional to u, is 
saved at 10 s intervals. This reporting interval is conserva- 
tive, as laboratory thermocouples can be sampled more than 
once pre second. In doing so, a time series of the "true" 
synthetic thermosyphon state is recorded to be used in a 
forecasting scheme. 



2.2 Forecast model 

The Ehrhard-Miiller (EM) system is a three-variable 
ODE deriv ed specifically to model b ulk flow in the ther- 
mosyphon (jEhrhard and Mulleij . [l990l : also see Appendix SI 
in the Supporting Information for an alternative derivation). 
Written in dimensionless form, the governing equations are 



dxi 
IF 

dX2 



a{x2 — xi) 



dt 

dX3 

IF 



- = 13x1 — X2{1 + Kh{\xi\)) — xiX'i 



X1X2 - x-i{l + Kh{\xi\)) . 



(2a) 
(2b) 
(2c) 



The state variable xi is proportional to the mass flow rate or 
mean fluid velocity, X2 to the temperature difference across 
the convection cell (ATs-g, measured between 3 o'clock and 
9 o'clock), and 2:3 to the deviation of the vertical tempera- 
ture profile from the value it takes during conduction; specif- 
ically, X3 oc {^ATw — AT6-12), where Ar6-i2 is the temper- 
ature difference measured between 6 o'clock and 12 o'clock. 
The parameter a is comparable to the Prandtl number, the 
ratio of momentum diffusivity and thermal diffusivity. Simi- 
lar to the Rayleigh number, the heating parameter /3 oc ATm 
determines the onset of convection as well as the transition 
to the chaotic regime. Finally, K determines the magnitude 
of variation of the wall heat transfer coefficient with veloc- 
ity. The functional form of that variation is determined by 
/i : R+ R+, where 



h{x) = 



55 3 I 20 4 



when X < 1 
when a; > 1 



(3) 



The interested reader is referred to Appendix SI in the 
Supporting Information for an explanation of this piece- 
wise form, which dif f ers sli ghtly from the original model of 
lEhrhard and Miillerl (|l990l ). 

Note th at when K = 0, the system is analo- 
gous to the iLorenj (|l963l ) system with geometric fac- 
tor (Lorenz's b) equal to one. The lack of a geometric 
factor in the EM system is due to the circular geom- 
etry of the convection cell. Lorenz equations have been 
widely used in nonlinear dynamics to stu dy chaos and in 
NWP as a model system for testing DA (Miller a nd Ghil . 
I994I: lYuen and Baul. 1 19991: I Annan and Hargrcavcsl. 12004: 



Evans et al] , 12004 lYang et al.1 . l2006l : iKalnav et al.l ■ l2007^ . 



When in the chaotic parameter regime, the EM system 
exhibits growing oscillations in the xi and X2 state vari- 
ables around their convecting equilibrium values until fiow 
reversal. In this system, the CCW rotational state is charac- 
terized by a;i > and X2 > 0, and the CW rotational state 
by xi < and X2 < 0. However, one should note that near 
a flow reversal xi and X2 can have opposite signs, because 
zero-crossings of the xi variable typically lag behind those 

of X2. 




Figure 2. An illustration of the basic predict-observe-update DA 
cycle. The EM model states (background forecast and analysis) 
are transformed into observations of the mass flow rate q by the 
observation operator (Eqn. (|4]l) for comparison with the truth. 
Here, using the 3D-Var algorithm and an assimilation window of 
5 min, the filter has satisfactory overall performance (scaled error 
~ 35%). Note the error spike around 135 min when the forecast 
and truth end up in different rotational states. The largest errors 
tend to occur at or near flow reversals due to inherent sensitivity 
near that transition and to the qualitatively different behavior of 
the different flow directions. 



The parameters found to match the simulated ther- 
mosyphon were a = 7.99, P = 27.3, and K = 0.148. The 
characteristic time and mass flow rate scales, used to trans- 
form the dimensionless model variables t' and xi into di- 
mensional time and "observations" of mass flow rate, were 
631.6 s and 0.0136 kg/s, respectively. The q scale is the one 
nonzero entry in the observation operator H, Eqn. (|4]). The 
above parameters were found using a multiple shooting algo- 
rithm explained in Appendix SI. 2 in the Supporting Infor- 
mation. Numerical integration of this autonomous ODE was 
performed with a fourth-order Runge-Ku tta method an d 
timestep 0.01 (corresponding to 6.316 s) in iMatlabI l|2009t ). 

2.3 Data assimilation 

DA is the process by which observations of a dynamical 
system are combined with forecasts from a model to estimate 
error covariances and calculate an optimal estimate for the 
current state of the system, called the analysis. The inherent 
difficulties are compounded by the fact that the forecaster 
uses an inexact forecasting model and never knows the true 
state of the dynamical system. The number of state vari- 
ables in a NWP model is typically 0(10'^) times larger than 
the number of observations. Nevertheless, the analysis be- 
comes the IC for a new forecast. The time interval between 
successive applications of the DA algorithm, i.e. the time 
between analysis steps (usually determined by the availabil- 
ity of observations but here allowed to vary), is called the 
assimilation window. The process is illustrated in Fig. (2] 

A variety of filters are capable of solving the DA 
problern . The canonical example is the Kalman filter (KF; 
lKalmanl.[l960l '). the optimal state estimation algorithm for a 
linear system. One of DA's first applications was to tra- 
j ectory estiniation and correction of missiles and rockets 
jSavelv et all . 1 19721 ). A number of nonlinear DA schemes are 
implemented in this study. In 3D variational DA (3D-Var; 
here 3D refers to the spatial dimensions for weather mod- 
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els), the background error covariance is estimated a single 
time, offline, prior to the data assimilation procedure. In the 
extended Kalman filter (EKF), background error is evolved 
according to the linear tangent model, which approximates 
the evolution of small perturbations about the trajectory. 
Ensemble Kalman filters (collectively EnKFs) use ensem- 
bles of forecasts to estimate the background error and bet- 
ter capture nonlinear behavior. The methods examined in 
this study were 3D-Var, the EKF, the ensemble square root 
Kalman filter (EnSRF) , and the ensemble transform Kalman 
filter (ETKF). Detailed descriptions of each method are in- 
cluded in Appendix S2 in the Supporting Information. A full 
review of DA is beyond the scope of the present pa per; for 
a com prehensive treatment, we refer the reader to iKalnavl 



3 Data assimilation experiments 

3. 1 Methods 

A perfect model experiment, in which the Lorenz equa- 
tions were used to forecast a synthetic truth created by 
the exact same system, was tested first but not included 
he re. We found ana lysis errors similar to t hose reported 
bvlYang et al.1 (|2006l ') (3D-Var and EKF) and lKalnav et al.1 
l|2007^ (ETKF), using the same model and tuning param- 
eters. This ensured that the DA algorithms were working 
before applying them to the synthetic thermosyphon data. 

As stated in Sec. 12.11 forecasts of the thermosyphon 
are made observing one scalar variable, the mass flow rate 
q oc x\. Gaussian noise with standard deviation equal to 
6 X 10^* kg/s, approximately 0.8% of the mass flow rate cli- 
matological mean, a/ {q^) = 0.075812 kg/s, is added to the 
synthetic truth to create observations. The relative mag- 
nitude of this error is comparable to that of experimental 
measurements. 

The EM model is used in the forecast step to integrate 
the analysis forward in time and create the new background 
forecast. The end results of applying DA are a background 
and analysis timeseries of xi,X2,X3, informed by both the 
timeseries of thermosyphon mass flow rate and the EM 
model dynamics. 

In this realistic forecasting scenario, where only limited 
information about the true state is available, the observa- 
tions of state variable q provide the only validation. For this 
reason, we calculate the forecast errors in observation space. 
These are given as root mean square error (RMSE), where 
RMSE = {Sq'^). The residual at a specific assimilation cy- 
cle is given hy Sq — q — Hx''. Here, x*" is the background 
forecast made by the model, and H : R'^ — >■ R is the linear 
observation operator 

H = [0.0136,0,0] (4) 

in units of kg/s. All errors are then scaled by the cli- 

matology of q. Analysis error is a common metric for assess- 
ing DA performance in perfect model experiments. In this 
study, however, we assert that background error is prefer- 
able. Analysis error in observation space, which will be small 
even for large assimilation windows, is not an appropriate 
metric for assessing model performance since it can disagree 
substantially with the background error. For example, 3D- 
Var in one experiment with a 10 minute assimilation window 



yielded analysis and background scaled errors of 0.08 and 
0.86, respectively. The analysis error would seem to indicate 
that forecasting is doing a good job, but the background 
error shows that background forecasts are essentially mean- 
ingless. The filter, however, accounts for this and weights 
the observations heavily over the background forecasts when 
producing the analysis. Since we are concerned with fore- 
casting, background error is a more representative metric. 

When applying DA to nonlinear systems, some type of 
covariance inflation is performe d to prevent fil t er div ergence 
due to error underestimation. iKalnav et al.l (|200'i1 ) found 
that a Lorenz forecasting model with a slightly different 
forcing parameter required a 10-fold increase in the mul- 
tiplicative inflation factor when using a 3 member EnKF. 
Model error is more pronounced for our forecasts, since the 
EM model is a reduced approximation of the numerically 
simulated thermosyphon. We relied upon additive and mul- 
tiplicative background covariance inflation to capture model 
error. Additive inflation was particularly important for the 
stability of the EKF and EnKFs. Additive noise provides 
a different exploration of dynamically accessible regions of 
state space, and it would be interesting to explore why ad- 
ditive versus multiplicative is preferred in certain cases, al- 
though this is beyond the scope of this paper. The speciflcs of 
how inflation was performed and tuned, and the parameters 
used are given in Appendix S2 in the Supporting Informa- 
tion. 

All EM and DA parameter tuning was performed us- 
ing a separate mass flow rate time series than was used for 
validation. Each DA algorithm was allowed 500 cycles to 
spin-up, and its performance was measured over the follow- 
ing 2500 cycles. Ensemble size in each case was set to 10 
members. 

3.2 Results 

With proper tuning, all DA algorithms were capable of 
synchronizing the EM model to observations of mass flow 
rate alone. As the assimilation window increased, scaled 
background error increased in a sigmoidal fashion, as ex- 
pected (see Fig. For assimilation windows up to 2.5 min, 
all DA algorithms have nearly indistinguishable errors. For 
assimilation windows between 3 and 6 min, 3D-Var performs 
noticeably worse than the other methods which remain in- 
distinguishable. Then, with assimilation windows greater 
than 6 min, the ensemble methods (EnSRF and ETKF) 
outperform the EKF noticeably. This is perhaps surpris- 
ing, at flrst glance, because the ensemble size is signiflcantly 
smaller than the dimension of the simulated thermosyphon 
state space (0(10®) variables). However, we know the ther- 
mosyphon dynamics effectively take place on the EM equa- 
tions' attractor (a manifold in three dimensions). The supe- 
rior performance of EnKFs here is likely due to the ensem- 
ble methods capturing nonlinear effects which dominate at 
larger windows. 

Following the historical SI score convention, scaled er- 
ror above 70% is considered a "usel ess" for e cast, while un- 
der 20% the forecast is "perfect" (|Kalnavl. 120021 '). Perfect 
forecasts for 3D-Var were found up to a 4 minute assimila- 
tion window, while the other methods (EnSRF, ETKF, and 
EKF) produced perfect forecasts with assimilation windows 
1.5 minutes longer. 
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2 4 6 8 10 

assimilation window (min) 

Figure 3. Background RMSE scaled by y' (q^) over 2500 as- 
similation cycles plotted for different DA algorithms and varying 
assimilation windows. As the window becomes larger, the error 
increases towards saturation. The lower dashed line in the main 
figure shows the limit of a "perfect" forecast while the upper de- 
marcates a "useless" forecast. 



A persistent spike in background error for the 5 minute 
assimilation window (Fig. |3]) is possibly due to that time 
being approximately the same as the characteristic period 
of oscillations in q (evident in Fig. [2| . We conjecture that 
this may lead to a type of resonance in the DA-coupled EM 
system which degrades DA performance. 

Besides these results pertaining to forecast skill, we also 
found that the DA algorithms infer thermosyphon dynam- 
ics which are absent from the EM model. In Fig. |4] we see 
the simulated thermosyphon's attractor o btained by b oth a 
time-delay embedding (Fig. |4ja) ; Alligoo d et al.l . Il996h and 
a projection of the EM analysis states to the X1-X2 plane 
(Fig.[4jb)). If the thermosyphon fluid flow stalls in the midst 
of a reversal, fluid in the bottom can quickly heat up while 
that in the top is cooled, leading to an unstable, strong 
temperature inversion. This causes the fluid to move very 
quickly in the reversed direction, but this new direction also 
ends up being unstable, and a new flow reversal can occur 
immediately. Absent DA, the EM model system does not 
exhibit this behavior of stalling followed by large swings of 
the trajectory. 

In the time-delay embedding (Fig. [H^a)), this phe- 
nomenon is exhibited by small loops in the trajectory as 
it moves near the convective fixed points. The flow stalls 
when the system state is near the conductive fixed point at 
the origin, then it swings wildly which brings it near the con- 
vective fixed point, but in such a way that it does not end 
up spiraling outward in the usual fashion as during a normal 
fiow reversal, as exhibited by the Lorenz equations. Instead, 
it quickly reverses again, which we call non-Lorenz behav- 
ior. This non-Lorenz behavior is further elaborated upon 
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Figure 4. Two views of the numerically simulated thermosyphon 
attractor. A time-delay reconstruction, using the monitored mass 
flow rate, is shown in (a). In (b), plotted points show xi and X2 of 
the EM analysis generated by EKF with an assimilation window 
of 120 s. Each is colored by the scaled background forecast error 
at that point. The delay time of 60 s used for (a) was chosen 
because it yielded a good approximation of the attractor in (b) . 
Note how in (a) trajectories that move through the far edge of 
either lobe create distinctive loops near the center of the opposite 
lobe. This is an example of dynamics which arc not present in 
the EM model without DA. It may explain the higher error for 
points in (b) at the far edge of each attractor lobe. See text for 
further description and Fig. [5] for another example. 



in Sec. 14.51 Forecast skill is worst at the far edges of the 
assimilated attractor (Fig. IHb)). This could be due to the 
wild swings of the EM trajectory after being ejected from 
the region of state space near the conducting equilibrium, 
or to the nonlinear dynami c al ins tabili ties at the edge of t he 
attractor found bv lPalmeil (|l993l ) and lEvans et all (|2004l ). 
We also explicitly show one of these stalled fiow rever- 
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Figure 5. The assimilated trajectory during a very large- 
oscillation, non-Lorcnz flow reversal. The EKF algorithm with 
a 30 s assimilation window was used. Color indicates the scaled 
background error. The state starts (a) in the CCW rotational re- 
gion and stalls near the conducting state for an extended period, 
causing fluid in the bottom to heat up, manifesting in an X3 that 
creeps towards with xi,X2 ~ 0, before the state swings quickly 
(b) through one oscillation in the CW rotational state. This is 
followed by one oscillation in the CCW state (c) before another 
stall near the conducting point and subsequent swing (b again) 
before settling into Lorenz-like oscillations (d). Note that the fil- 
ter only observes the xi variable but is able to reconstruct the 
dynamics in the full state space. 

sals in Fig.O where we plot the EKF-assimilated EM trajec- 
tory using a 30 s assimilation window. When the fluid stalls, 
the X3 variable moves closer to (i.e. Ar6_i2 increases) 
while xi and X2 (proportional to q and ATs-g, respectively) 
are approximately 0, reflecting the growing temperature in- 
version while the fluid remains stationary. When the fluid 
starts to move, the assimilated trajectory swings wildly to 
the left attractor lobe (CW rotational state), then right 
(CCW rotational state). The trajectory undergoes another 
stall-swing cycle before flnally resuming Lorenz behavior, 
where the trajectory spirals outward from the CCW convect- 
ing equilibria. This contrasts the Lorenz and EM model dy- 
namics, for which large deviations in the system state from 
convecting equilibrium are driven close to the other con- 
vecting equilibrium during a flow reversal, which stabilizes 
the system. See also Sec. 14.51 Fig. [T] and the accompany- 
ing discussion. This result remains unchanged for the other 
DA algorithms also using a 30 s assimilation window. The 
inference remains using EKF and a 60 s assimilation win- 
dow, but the trajectory appears much noisier, leading us to 
believe that this is due to the rapid update. With larger as- 
similation windows, the trajectory becomes uninterpretable 
as error in the unobserved variables increases. 



4 Flow reversal experiments 

4.1 Experimental setup 

For the purpose of flow reversal forecasts, we picked a 
single DA algorithm and assimilation window. In this Sec- 
tion, all analyses were generated by the extended Kalman 
filter and an assimilation window of 30 s. This interval corre- 
sponds to 5 time steps of the model and is shorter than that 



used in lYang et al] (|2006l ) and lKalnav et all (|2007l '). The fol- 
lowing could certainly be repeated using other algorithms, 
observations, and assimilation windows, but this was beyond 
the scope of this paper. The fiow reversal tests in Sec. 14.31 
and the residency time forecasts in Sec. 14.41 were tuned and 
validated on separate analysis timeseries. The length of the 
tuning and validation timeseries were approximately 39 and 
93 days, respectively. 



4-2 Occurance of flow reversals: traditional explanation 

The first explanation of the mechanism respo nsible for 
fiow rever sals was presented by IWelanderl (| 19671 ) and re- 
peated by ICreveling et al] (|l975h . Welander, who was also 
the first to discover that thermosyphons exhibit aperiodic 
oscillatory behavior, explained the instability of steady con- 
vecting fiow by considering a thermal anomaly or "warm 
pocket" of fiuid. For low heating rates, the convecting equi- 
librium is stable because viscous and thermal dissipation are 
in phase, thus an increase (decrease) in fiow rate leads to an 
increase (decrease) in friction and a decrease (increase) in 
buoyancy, and such perturbations are damped out. At higher 
heating rates, the warm pocket is amplified with each cycle 
through the loop due to out of phase viscous and thermal dis- 
sipation's. Welander explained that when the warm pocket 
emerges from the heating section and enters the cooling sec- 
tion, it feels a greater buoyant force than the surrounding 
fiuid and accelerates, exiting the cooling section quickly, giv- 
ing it less time to radiate away its energy. As the pocket 
moves into the section with warm boundary, the buoyant 
force it experiences is again higher than normal, so now the 
pocket decelerates and passes slowly through the heating 
section, gaining more energy. This positive feedback effect 
causes the pocket to grow hotter and larger with each pass 
through the loop. These oscillations in the fiuid temperature 
and velocity do not grow unhindered, however. The pocket 
eventually becomes large and hot enough that its descent 
towards the heating section is stopped entirely by its own 
buoyancy. Without movement, the pocket dissipates, but its 
remnant heat biases new rotation in the opposite direction, 
and the fiow reverses. 

In the Lorenz and EM systems, this feedback is embod- 
ied in the spiraling repulsion of trajectories from the unsta- 
ble convecting equilibria at the center of each lobe or wing of 
the attractor before moving to the other lobe. Because the 
growth of oscillations is an important component to the fiow 
reversal process in both the CFD simulated thermosyphon 
and EM system, we define here what is meant by oscilla- 
tion amplitude in each case. In the CFD simulated ther- 
mosyphon, it is the maximum distance of the system state 
from the nearest convecting equilibrium, where system state 
is understood to mean the state of the entire temperature 
and velocity fiow fields in the CFD simulation. When consid- 
ering the DA-generated EM analysis, the kth xi oscillation 
amplitude a;™"" is defined as the maximum amplitude 

xi^'^^ = ma.x\xi(t)\ (5) 
teT 

where T — [io*'i^/*'l is the time interval of the kth oscilla- 
tion. 
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4-3 Flow reversal forecasting methods 

Three separate tests were developed to predict, at each 
assimilation step, whether a flow reversal would occur within 
the next oscillation period (approximately 11 min), here 
taken to be within the next 20 DA cycles. See Sec. 14.61 and 
Appendix S3 in the Supporting Information for a description 
of how the tunable parameters were chosen. 



4-3.1 Lead forecast The simplest test forecasts a flow re- 
versal whenever the background forecast changes rotational 
state. Note that to forecast a flow reversal occurring in the 
future, the background forecast started from the most recent 
analysis IC provides our only information about the system's 
future state. Ignoring the three-dimensional nature of the 
state space, a flow reversal is forecast whenever xi crosses 
through zero. Note also that the forecast is unable to pre- 
dict flow reversals that occur beyond the lead time, and that 
lead forecast quality quickly degrades as the lead time is in- 
creased. We impose a limit on the number of assimilation 
cycles to look ahead, Aiead ~ 7, so that the algorithm does 
not trust forecasts too far in advance. 




— threshold 

"^"o 20 40 60 80 100 120 140 160 180 200 

time (min) 

Figure 6. Correlation test: Whenever the slope (green) of the 
[x2,xi]'^ correlation exceeds a threshold (red), a flow reversal is 
forecast. Correct positive predictions are shown as filled circles 
and false positives as open circles. The starred point corresponds 
to the inset, which shows how correlation is computed as the slope 
of the least squares fit (green line) of previous analysis points, and 
the arrow shows the direction of the trajectory. Note that each 
analysis cycle where the correlation fails to exceed the threshold 
counts as a correct negative forecast (not shown). There are no 
false negatives, i.e. misses, in this timescries. Here pcorr and Acorr 
are the same as for Tab. [T] 



4-3.2 Bred vectors An ensemble of perturbed states form- 
ing a small ball around the analysis can be used to represent 
uncertainty in the IC. A nonlinear system will dynamically 
stretch and shrink such a ball around its trajectory as i t 
moves through the attractor (jPanforth and Yorkd . l2006h . 
Small perturbations to points on a trajectory are integrated 
forward in time, and the differences between perturbed and 
unperturbed solutions are called bred vectors (BVs). Here, 
the rescaling amplitude is 0.001 and the integration time 
coincides with the 30 s assimilation window. 

The average BV growth ra t e is a useful measure of lo - 
cal instabilities (|HofFman et al.l . |2009| ). lEvans et all (|2004 ). 
studying perfect-model forecasting of the Lorenz system, set 
a BV growth rate threshold which accounted for 91.4% of the 
observed flow reversals (hit rate). Our BV test simply fore- 
casts a flow reversal whenever the average BV growth rate 
over the previous assimilation window exceeds a threshold, 
PBv = 0.6786. 



4-3-3 Correlation The final test uses the fact that flow 
reversals are suspected to be caused by out of phase viscous 
and thermal dissipation. Since the friction term grows with 
fluid velocity oc xi and the thermal dissipation grows with 
the size of the temperature anomaly, related to X2, we exam- 
ined the correlation between those two variables over a tun- 
able number of previous analysis cycles. Specifically, when 
the slope of the least-squares linear fit of Acorr ~ 18 previous 
analysis points [x2{t—i),xi{t—i)]'^ for i = 0, 1, . . . , (Acorr — 1) 
exceeds a threshold pcorr = 1.42, a flow reversal is forecast. 
See Fig.[6]for an illustration of this process. Interestingly, in- 
creasing autocorrelation of the state seems to be a universal 
property of many sy s tems in advance of critical t ransitions 
IScheflfer et all . l2009l : ICotilla-Sanchez et al.l . [20]l ). 



4-4 Forecasting residency times in the new rotational state 

We found that the analysis' xi oscillation amplitude 
preceding each flow reversal is correlated with the duration 
of the following rotational state, shown in Fig. [7] We refer 
to these durations between flow reversals as residency times. 
Residency times are observed at discernible "steps" corre- 
sponding to integer numbers of oscillations. This correlation 
makes the xi oscillation amplitude a plausible predictor for 
residency time in the new rotational state. 

Furthermore, the average BV growth rate measured 
over the assimilation window preceding that extremum fol- 
lows a clear gradient in Fig. [3 the growth rate increasing 
with oscillation amplitude. The BV growth rate gradient 
implies that more unstable system states precede longer 
residency times in the next rotational state. Outliers with 
^max > -j^^ g result in shorter residency times than expected 
from making similar plots to Fig. [7] for the pure Lorenz and 
EM systems (not shown). In the Lorenz and EM systems, 
the steps continue to move upwards with xi oscillation am- 
plitude. The discrepancy is due to the non-Lorenz behavior 
that was mentioned at the end of Sec. 13.21 

Our residency time prediction algorithm proceeds as fol- 
lows. When a flow reversal is forecast by one of the methods 
described in Sec. 14.31 the algorithm first calculates x™^" as 
defined by Eqn. ((Sjl for the presently occurring oscillation. 
The algorithm uses only the analysis and lead forecast data 
available at the time the flow reversal test is triggered when 
estimating x™'"'. The algorithm then examines the residency 
times of all flow reversals which followed an a;™^" in the in- 
terval (2:5"^" - 0.5, a;']"^" + 0.5). These ordered pairs of am- 
plitudes and residency times are drawn from the training 
timeseries. From the relative abundance of residency times 
in this sample, we assign a probability to the number of flow 
oscillations in the forthcoming rotational state. (See the in- 
set histogram in Fig.[71) The categories are restricted to 1-6 
oscillations (a duration of 7 oscillations, shown in Fig. [T] is 
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Figure 7. Residency time in the new rotational state is plotted versus the amplitude of xi (proportional to the mass flow rate) at the 
last extremum before flow reversal. This amplitude is calculated from the EM model EKF analysis of thermosyphon observations, using 
a 30 s assimilation window. This figure contains over 39 days of simulated flow and 1796 flow reversals. Points are colored by the average 
BV growth rate over the preceding assimilation window, showing a BV growth rate gradient that increases in the positive direction along 
both axes. Inset, upper left: A timeseries corresponding to a single point in the scatter plot, marked with a black cross. The xi oscillation 
amplitude and subsequent residency time are denoted by the vertical and horizontal arrows, respectively. Inset, upper right: A histogram 
showing the likelihood of residency times following an x^^^ G (11, 12). This is the interval that we would consider for an xi oscillation 
amplitude of 10.5 preceding flow reversal. The most likely residency time is about 23 min or 2 oscillations, the middle "step" for the 
histogram range. 



observed exactly once in the training timeseries, so it was 
considered too rare an event to merit a category). The typ- 
ical residency times corresponding to 1, 2, 3, 4, 5, and 6 
oscillations are taken to be 11.48, 23.09, 33.72, 44.38, 55.11, 
and 66.08 minutes, respectively; the oscillation category as- 
sociated with a given residency time is taken to be that 
with the closest time in this list. This algorithm generates a 
probabilistic forecast from the relative abundance of points 
in each oscillation category. An example output would be 
20%, 40%, 30%, and 10% chance of 1, 2, 3, and 4 oscilla- 
tions in the next rotational state and zero probability of 5 
or 6 oscillations. 



4-5 New details regarding the flow reversal mechanism 

Not all flow reversals occur when the system reaches the 
same flow oscillation amplitude, nor do all rotational states 
last the same amount of time. During a flow reversal, the 
fluid motion stalls after hot fluid extends across the entire 
heating section into the cooling section (see Sec. 14.21 and 
Fig-O- The magnitudes of this hot "tongue" and, likewise, 
the opposite cold tongue affect the stability of the system as 
it reverses. If the oscillation is small, it will mostly dissipate 



before the new rotational state is entered, bringing the tem- 
perature profile close to that of conduction. This is a highly 
unstable equilibrium, since the vertical temperature gradi- 
ent builds until the fluid in the bottom is much hotter than 
the fluid above (illustrated in the analysis in Fig. [5]). When 
the fluid begins to rotate, it accelerates rapidly. The large 
amount of heat carried by the fluid brings the system state 
far from the convecting equilibrium. If the oscillation is large 
(corresponding to a large deviation from convecting equilib- 
rium in temperature and velocity), remnant warm and cool 
areas will be present in the top and bottom sections of the 
loop, respectively. These stabilize the new rotational state 
near its convective equilibrium. The resulting duration is 
longer since the instability requires more time to grow be- 
fore causing the next reversal. These two situations are illus- 
trated in Fig.|S]and explain the trend in the Lorenz region of 
Fig. [7] Animations of the simulated temperature fleld during 
flow reversal are consistent with this explanation 0. 

We believe that the behavior in the extremely large os- 
cillation, non-Lorenz region, where x™^" ^ 14.5 and shown 

^ A movie similar to the case shown in Fig.|8]is provided online at 
|http : //www.uvm. edu/-kharris/therinosyphon/T-Ra-18000-new.mp4 
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Figure 8. Temperature profiles before and after two fiow reversals for the chaotic regime in units of Kelvin. For ease of visualization, 
the radius ratio is reduced to 3. Ra = 1.8 X 10*. The amplitude of the oscillation shown in (a) is relatively weak; the temperature 
profile quickly approaches (b) that of conduction (c) where it heats up significantly before the reversal. The extreme instability of the 
conducting state, near time 20 min, produces a large oscillation in the CW direction (d), immediately causing another flow reversal back 
to CCW. As the system enters the new rotational state, remnant heat stabilizes the flow [contrast (e-f) with (b-c)], necessitating a longer 
residency in the new rotational state while the instability grows (g-i) (note that for the radius ratio of 24 no more than 7 oscillations are 
ever observed). 
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in Fig. [T] is caused by excessive remnant thermal energy 
after flow reversal. Although the temperature distribution 
present after a flow reversal is configured in a way that sta- 
bilizes the flow in the new rotational state, the very large 
magnitude of the temperature field is a competing, destabi- 
lizing factor that dominates as xT'^^ increases into the non- 
Lorenz region. This leads to shorter durations in the new 
rotational state before a second flow reversal occurs. 



(a) Lead forecast, n=175592 



Observed 



Yes 



No 



Feast 



Yes 4472 



744 



No 



13 



170363 



(b) BV test, 11=121120 



4-6 Flow reversal forecast skill 

The results of the three tests are presented in Tab. [1] 
as two-by-two contingency tables. Shown in Tab. [2] are the 
threat score (TS), f alse alarm ra tio (FAR), and probability 
of detection (POD) (| Wilksl. [l995l ) . Given a non-probabilistic 
yes/no forecast with a hits, b false alarms, c misses, and d 
correct negatives for a total of n events, these are defined as 
TS=a/{a + b + c), FAR=b/{a + b), POD^a/{a + c). Because 
flow reversals are relatively rare events, the hit rate {a + d)/ri 
would be dominated by correct negatives. Instead, TS is 
chosen as an appropriate overall performance metric since 
it disregards these frequent negative events and takes into 
account both false alarms and misses. 

There are trade-offs among the various skill scores 
for each flow reversal test. Tuning the reversal tests then 
amounts to multiobjective optimization, attempting to max- 
imize TS, RPS-avg, and RPS-med (the skill scores used for 
residency time forecasts, defined in Sec. l4.7|) . minimize FAR, 
and maintain POD above 95%. The goal was to tune each 
method to all-around good performance, for both reversal 
occurrence and residency time forecasts. To guide the pro- 
cess, plots of the skill scores were made for different tuning 
parameters, but the final tuning was performed ad hoc. In 
Appendix S3 in the Supporting Information, Fig. S3-2 shows 
one of these tuning experiments with the final parameters 
chosen appearing in the center of each subfigure. 

Considering TS alone, the lead forecast performed best, 
followed by the correlation test, with the BV test performing 
poorest. The BV test also had a very hi gh FAR, leading u s 
to conclude, in contrast to the results of lEvans et al.l (|2004l ') 
for a perfect model experiment, that BV growth rate is a 
poor overall predictor of flow reversals in a realistic ther- 
mosyphon. On the other hand, the correlation test had the 
lowest FAR while maintaining a high TS, but this comes 
at the price of more misses, resulting in a lower POD. The 
reasonable performance of the correlation test in all areas 
lends circumstantial evidence to the claim that out of phase 
dissipations are indeed the cause of flow reversals. 

The flow reversal occurrence tests are triggered in dif- 
ferent situations, leading to variation in how far in advance 
flow reversals are detected, the "warning time". Warning 
times were only computed for hits, i.e. forecast flow reversals 
which were observed to occur. The lead, BV, and correlation 
tests had average warning times of 175, 217, and 304 s re- 
spectively. Histograms of these warning times are presented 
in Appendix S3 in the Supporting Information. 



4-7 Residency time forecast skill 

A naive way of forecasting residency times would assign 
each possible outcome a probability equal to that measured 



Observed 



Yes 



No 



Feast 



Yes 4383 



3203 



No 



102 



121258 



(c) Correlation test, n=174925 



Observed 



Yes 



No 



Feast 



Yes 3540 



239 



No 



945 



170201 



Table 1. Contingeney tables for the three flow reversal tests. The 
parameters used were: pcorr = 1.42, Acorr = 18, PBV = 0.6786, 
and Aioad = 7. For a detailed description of the tuning, see the 
text and Fig. S3-2 in Appendix S3 in the Supporting Information. 



Method 


TS 


FAR 


POD 


RPS-avg 


RPS-med 


lead 


86 


14 


>99 


71 


87 


bred vector 


57 


42 


98 


67 


86 


correlation 


75 


6 


79 


58 


74 



Table 2. Skill scores for flow reversal categorical forecasts (TS, 
FAR, POD) and residency time probabilistic forecasts (RPS-avg, 
RPS-med) for the three tests. Note that for TS and POD a perfect 
score is 100%, while a perfect FAR is 0%; RPS scores should be 
interpreted as a percent improvement over climatology, so that 
any RPS > is an improvement. 



from the climatology. In our case, this would amount to us- 
ing the marginal distribution of oscillation occurance. How- 
ever, our method also takes into account the a;™^" before 
the flow reversal (i.e. the joint distribution of events by os- 
cillation occurance and xf^'^^) , which we have shown contains 
important information about the number of oscillations that 
the system will undergo in the new rotational state. So, we 
compare our method to climatology using a ranked probabil- 
ity skill score (RPS, see [Wilks . 19951 ). This is only computed 
in the case of hits. We actually computed two variants, by 
taking either the mean (RPS-avg) or median (RPS-med) of 
the ranked probability scores for each reversal event when 
computing the skill. The results are illustrated in Tab.[51 The 
lead forecast test performs best, followed by the BV test and 
the correlation test. Unsurprisingly, the flow reversal tests 
with smaller warning times performed better when making 
residency time forecasts. Because there is more information 
about the system state immediately preceding a flow rever- 
sal if the warning time is small, the residency time forecast is 
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better informed. The magnitude of the improvement over cli- 
matology is large for all methods. The RPS scores in Tab. [2] 
are simila r to or better than those for probabilistic forecasts 
in NW P (jTippett and Barnstonl . l2008l : [Poblas-Reves et al.l . 
|2000D . 



5 Conclusion 

DA was shown to be an effective way of coupling a sim- 
plified model to CFD simulations of the thermosyphon. Al- 
though background forecast errors were always larger than 
observational noise, climatically scaled background error 
was small for reasonable assimilation windows. Proper tun- 
ing of multiplicative and additive inflation factors was essen- 
tial for avoiding fllter divergence and achieving low forecast 
error. All of the DA methods used in this study accurately 
capture the behavior of the thermosyphon with short as- 
similation windows. Among the DA methods, the ensem- 
ble methods show advantages over 3D-Var and EKF with 
longer assimilation windows, when nonlinear error growth 
becomes important. With frequent analysis update, DA can 
reveal non-Lorenz behavior in the thermosyphon even with 
the EM (Lorenz-like) model. 

Three different predictors of flow reversals were pro- 
posed and tested wi th reasonabl e succ ess. In comparison 
with the two rules in lEvans et all \2004 ] for predicting the 
behavior of the Lorenz trajectory, the BY growth rate is a 
useful measure of the EM model's dynamical instabilities, 
but it does not perform well on its own as a predictor of 
flow reversals. Finally, the amplitude of the final oscillation 
in the current rotational state was found to be correlated 
with the residency time in the following rotational state, 
and we provide a physical explanation for this phenomenon, 
elaborating on the details of fiow reversals. Oscillation am- 
plitudes were then used to create probabilistic forecasts of 
those residency times with significant improvements over cli- 
matology. 

A laboratory thermosyphon device is in construction. 
The next stage of this research will apply similar meth- 
ods to forecasting the system state, flow reversals, and resi- 
dency times using 3D numerical flow simulations. Spatially- 
aware DA techniqu es, such as the Local Ensemble Trans - 
form Kalman Filter (jKalnav et al.1 . l2007l : [Hunt et al.l , l2007^ ■ 
could be applied to finite-volume or finite-element mod- 
els to study the spatial structure of the fiuid flow and er- 
ror growth. These imperfect model experiments could be 
used to compare t he relative perform ance of other DA al- 
gorithms (4D-Var. lKalnav et al.l. | 2007l). svnchroni zation ap- 
proaches (adaptive nudging, see lYang et al.1. |2006|) . and em- 
pirical correction techniqu es (|Danforth et al.l . 2007l : lLret al.l . 
l2009l : lAllgaier et all . |2012| ) . 

Although the thermosyphon is far from representing 
anything as complex and vast as Earth's weather and cli- 
mate, there are characteristics our toy climate shares with 
global atmospheric models. Sophisticated atmospheric mod- 
els are, at best, only an approximate representation of the 
numerous processes that govern the Earth's climate. Global 
weather models and the EM model both parameterize flne- 
scale processes that interact nonlinearly to determine large- 
scale behavior. Clouds and precipitation are sub-grid scale 
processes in a global weather model, and the correlations 



for the heat transfer and friction coefficients are parame- 
terizations of fiuid behavior on a finer scale than can be 
dealt with in the reduced model. Cloud formation is only 
partly understood, and moist convection is an area of ac- 
tive research where som e models bear similarit ies to the EM 
model considered here (jWeidauer et al.l . l201lf ) . 

The methods we use to forecast the toy model are also 
similar to the methods used for global geophysical systems. 
Both require state estimation to find the IC from which to 
generate forecasts. Also, when forecasts are made in either 
system, climatology and dynamically accessible regimes are 
often more important than specific behavior: the occurrence 
of fiow reversals for the thermosyphon; periodic behavior 
such as the El Niiio Southern Oscillation, and statistics such 
as globally and regionally- averaged temperatures and their 
effects on rainfall, ice cover, etc. for climate. Each of these is 
a statistic that must be post-processed from the model out- 
put. To meet these global challenges, many tools are needed 
in the modeling toolbox. These techniques may also be useful 
for forecasting sociotechnological systems which are rapidly 
gaining importance as drivers of human behavior. In this 
way, toy models can provide us with insights that are appli- 
cable to the important scientific problems of today. 
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APPENDIX SI: Ehrhard-Muller model 

Sl.l Derivation 

Following the d e rivatio ns of I Gorman et all (|l986l ) and 
lEhrhard and Mulleij (|l99d ). we consider the forces actins 
upon a control volume of incompressible fluid in the loop. All 
fluid properties are cross-sectionally averaged, and the radial 
components of velocity and heat conduction within the fluid 
are neglected. The fluid velocity u — u{t) is assumed to be 
constant at all points. Applying Newton's second law, the 
sum of all forces on the control volume must equal its change 
in momentum: 

Fp + Ff + Fg = pnr'^Rdfl,^ (Sl-la) 
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where 



Ofp 

Ff = —pnr^Rd(j)f^ 

Fg = — pnr^ Rd(j} g sin (j) . 



(Sl-lb) 

(Sl-lc) 
(Sl-ld) 



The angular coordinate and loop dimensions r and R are 
defined in Fig. 1 in the main text, g is the acceleration of 
gravity, p is the fluid density, u is velocity, and p is pres- 
sure. The total force in Eqn. (|Sl-la|) is comprised of the net 
pressure (Fp), friction from shear within the fluid (Ff), and 
the force of gravity {Fg). The pressure term, Eqn. (|Sl-lb|) . 
is the volume times the pressure gradient. The friction term, 
Eqn. (|Sl-lc|) . is written in this form in order to simplify the 
analysis; all frictional effects are contained in which will 
depend on fluid velocity, to be discussed later. 

Before we write the momentum equation, it is conve- 
nient to apply the Boussinesq approximation, which assumes 
that variations in fluid density are linear with temperature. 
In other words, p — p{T) « po{l — ^{T — To)) where po is the 
reference density, 7 is the coefficient of volumetric thermal 
expansion, and To = \{Th + Tc) is the reference temper- 
ature. The Boussinesq approximation also states that the 
density variation is insignificant except in terms multiplied 
by g. Thus, the density p is replaced by po in all terms of 
Eqn. (|S1-1|) except gravity, Eqn. HSl-ld|l . Using the Boussi- 
nesq approximation, gathering terms, and dividing out com- 
mon factors gives the momentum equation 

(Sl-2) 

Integrating about the loop, the momentum equation is sim- 
plified because u and are independent of 4> and other 
terms drop out due to periodicity. 



po 



du 
It 



po7g 

2ti 



i T sin — 



(Sl-3) 



We now must account for the transfer of energy within 
the fiuid, and between the fluid and the wall. All modes 
of heat transfer are neglected except c onvection, which 
is a valid approxim a tion w hen r <^ R (jWelandei) . 1 19671 : 
lEhrhard and Miilled . Il99d ). The energy rate of change 
{D/Dt is the material derivative with respect to time) in 
the control volume is 



p,^r Rd^Cp— = po.r Rd4>Cp{— + - — 



(Sl-4) 



which must be equal to the heat transfer through the wall 



AQ = ~TTrR dct> {T - 



(Sl-5) 



where Cp is the speciflc heat of the fluid, is the heat 
transfer coefficient, which depends on velocity, and r„ is the 
temperature at the wall. Combining Eqns. (|S1-4|) and (|S1-5|) 
gives the energy equation 



dT udT 
'dt ^ R'd^ 



K 

paCp 



-{T-T^ 



(Sl-6) 



Together, Eqns. HS1-3|I and (|S1-6|) represent a simple model 
of the flow in the loop. 

The transport coefficients and characterize the 



interaction between the fluid and the wall. They are defln ed 
by the constitutive relations (|Ehrhard and Miilleij , Il990l ) 



hw = hu,o [1 + Kh(\xi\)] 
1 



(Sl-7) 

fw = ^pofwoU, (Sl-8) 

where a;i oc u is the dimensionless velocity. The function 



in Eqn. (|Sl-7fl determines the velocity dependence of the 
heat tr ansfer coefficient, which va ries as u^''^ for moder- 
ate u (|Ehrhard and Miilled . Iiggd ). We introduce the fit- 
ting polynomial p{x) = 44/9 - 55/9 a;^ -I- 20/9 x"* to en- 
sure that /im is analytic at xi = 0. This piecewise defini- 
tion causes h{x) to vary as p{x) for a; 1 and x^''^ for 
3; > 1. Eqn. (tSTSll gives the frictional deceleration of the 
fluid when \u\ > 0, and the po/2 term is retained to sim- 
plify the final solution. Dimensionally, /„, is an acceleration 
(m/s^) and is power per unit volume per unit temper- 
ature (W/m'^K). These coefficients fe mn 1 fwn . and K must 
be estimated from experiments (e.g., lEhrhard and Miilleij . 
I1990I : IWelandeJ. 1 19671 : iGorman et all . 1 19861 ) or from other 
empirical means. In Sec. IS1.2I we describe the empirical 
met hods used for parame t er est imation. 

lEhrhard and Miilleij (|l990l ) solved the system of two 
coupled, partial differential equations (Eqns. (|S1-3|I and 
(|S1-6|) ) by introducing an infinite Fourier series for T. The 
essential dynamics can be captured by the lowest modes, 
i.e., 

T{(t>,t) = Co{t) + S{t) sin (f) + C{t) cos (t>. (Sl-10) 

Because this form of T separates the variables </> and t, the 
problem is transformed into a set of ordinary differential 
equations. Substituting Eqn. (|S1-10|) into Eqn. ()S1-3|) and 
integrating gives the equation of motion for it. Similarly, 
Eqn. (|S1-6[I is integrated by f d<f) sin cj) and § d(f) cos to sep- 
arate the two temperature modes 5* and C. The system is 
written in dimensionless form 



dxi , 

^^I3x,- 
dt' ^ 

dx-j 
dt' 



— X1X2 



-xi) (Sl-lla) 
X2{l+Kh{\xi\))-XlX3 (Sl-Ub) 
-X3{l+Kh{\xi\)) (Sl-llc) 



where the following linear transformations have been made 
to create dimensionless variables 



t = 



Xi = 



poCp 
poCp 

Rhw(\ 



1 

3:2 = 7: 



PoCp-yg 



2 RhwQ fwQ 



ATs 



X3 



1 Pocpig 

2 RhyjQ fwQ \ 



(Sl-12) 



Physically, xi is proportional to the mean fiuid velocity, 
X2 to the temperature difference across the convection cell 
or AT3_9 (between 3 o'clock and 9 o'clock), and X3 is pro- 
portional to the deviation of the vertical temperature pro- 
file (characterized by the temperature difference between 6 
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o'clock and 12 o'clock, AT6-12) from the value it takes dur- 
ing conduction. 

The parameter a = \poCpfwQ /h^^ is comparable to 
the Prandtl number, the ratio of momentum diffusivity and 
thermal diffusivity. Similar to the Rayleigh number, the 
heating parameter 

/3=-^^^AT„ (Sl-13) 

determines the onset of convection as well as the transition 
to chaotic behavior. 

Although the previous derivation assumes a 3D geome- 
try, the CFD simulations described in Sec. 2.1 of the main 
text, were performed in 2D. A 2D geometry corresponds 
to infinite concentric cylinders as opposed to the quasi-lD 
torus. Due to cross-sectional averaging, the EM equations of 
motion (|S1-11|) are the same in 2D or 3D; the change may 
be realized by letting vrr^ — >■ 2r in Eqns. (|S1-1|) . (|S1-4|) . and 
(|S1-5|) and carrying out the rest of the derivation. The only 
differences arise in the non-dimensional transformations and 
parameters, which were empirically determined by a multi- 
ple shooting algorithm (see Sec. IS1.2|) . 



0.2 




% 100 200 300 400 500 600 

time (min) 



Figure Sl-1. Results of the multiple shooting algorithm. The 
entire training timeseries is shown. Starred points are the CFD 
thermosyphon mass flow rate data used for validation. Also shown 
are the trajectories of model integrations over each shooting win- 
dow, which are indicated by the dashed vertical lines. In one 
shooting window, near 515 min, the model does not match the 
data because the optimization has not found a good IC (it is near 
the highly unstable conducting equilibrium) even though the pa- 
rameters are acceptable. 



SI. 2 Parameter estimation 

Before any forecasting, the parameters matching the 
EM model to the thermosyphon simul ation needed to be 
determined. lEhrhard and Miillerl (| 19901 ) used experimental 
measurements to determine the correlation coefficients for 
friction, /i^q, and heat transfer, hw^ and K. They achieved 
this by opening the loop at </> = 7r/2 and providing a devel- 
oped flow with adjustable velocity. By measuring the pres- 
sure loss (oc fw) and heat transfer across the loop for a 
range of velocities, they were able to find the correlation 
coefficients using regression. We were unable to accomplish 
this with a CFD simulation of an open-loop geometry. 

Instead, parameter estimation was formulated as a mul- 
tiple shooting problem. Shooting methods minimize the er- 
ror in an ODE trajectory relative to data by optimizing over 
all possible initial conditions and parameter space. Multiple 
shooting is a shoo ting method suitable for chaotic ODEs 
jBaake et al.1 . Il992l ). It overcomes the sensitive dependence 
on initial conditions by partitioning the data set and solv- 
ing the shooting problem on those subsets of the data, aug- 
mented by continuity conditio ns. We used the nonlinear least 
square optimizer Isqnonlin in iMatlabI (|2009l ') to perform the 
minimization and relaxed the continuity constraints. The 
model parameters which were tuned were a, P, and K. How- 
ever, we also needed a way to determine the time and ve- 
locity scales to convert the dimensionless variables t' and xi 
to their observed, dimensional values t and q. These scales 
change as the other parameters are varied, so these were 
incorporated into the variables of the optimization. 

The results of the multiple shooting algorithm are 
shown in Fig. lSTTI and Table [FTT] 

APPENDIX S2: Data assimilation algorithms 

S2.1 Kalman Filter (KF) 

The KF is well-known and widely used in linear DA 
and control problems. Although the thermosyphon is highly 



parameter value 

a 7.99 

/3 27.3 

K 0.148 

t scale (s) 631.6 

q scale (kg/s) 0.0136 

Table Sl-1. Final parameters used in the data assimilation 
scheme. The first three are the dimensionless parameters of the 
model, and the final two are used to rescale the dimensionless 
time and mass flow rate. 

nonlinear, the linear update equations are similar to those 
of the nonlinear algorithms used for this experiment. The 
KF attempts to assimilate observations and forecasts for a 
process of the form 

= Wxl._i . (S2-1) 

In this case, is the true state, which advances in time 
according to the linear process W, which is unknown but 
approximated by the model M. Subscripts index the time 
step. Using the model, the analysis from the previous time 
step is integrated to generate the background forecast for 
the current time step 

xt = Mx^i (S2-2) 

where M is the linear model, x" is the old analysis, and x*" is 
the background. Because M is only an approximation of W, 
a perfect initial condition will not lead to a perfect forecast, 
so 

X* = MxLi + ^1 (S2-3) 

where the model errors have covariance Q (usually as- 
sumed to be constant in time) and are written on the right 
hand side for convenience. When deemed unnecessary, time 
subscripts are left out. 

Given an observation y and background forecast x*", 
the KF finds the optimal way to combine them into the 
analysis x", the best guess of the current state. This becomes 
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the IC when forecasting with the model, Eqn. (|S2-2[) . In an 
operational context, we usually cannot observe every state 
variable. If x e R'^ and y G k'*^ then M < N {in NWP 
M <C N), so we define the observation operator H : — > 
R*^ that takes the background forecast from the model state 
space into the observation space. This serves two purposes: 
first, it avoids extrapolation of observations to gridpoints 
in state space; and second, it enables us to interpret our 
forecasts by comparing them directly to observations. For 
the thermosyphon, H is linear, so we write it as H, but this 
is usually not the case for the observations in NWP, e.g., 
satellite radiances and radar refiectivities. 

The complete application of the KF consists of a fore- 
cast step 

xt = Mx^i (S2-4a) 
Bk = MAfe-iM^ + Q (S2-4b) 
and an analysis step 

x^ = x^+Kfc(yfc-Hx^.) (S2-5a) 
Afe = (1 - KfcH)Bfc (S2-5b) 

with the Kalman gain given by 

Kfc =BfcH^(HBfcH^ + R)-\ (S2-6) 

The forecast equations create the background forecast 
and update the background error covariance. The new back- 
ground error covariance is the old analysis error integrated 
forward plus the model error Q. In the analysis step, this 
background forecast is incremented by the gain times the in- 
novation (y — Hx*") to produce the analysis. The difference 
between the analysis and the background is referred to as 
the analysis increment; statistical properti es of these incre- 
ments can be used to reduce mode l error jPanforth et al.l . 
I2OO7I : iDanforth and Kalnavl. l2008bH ah . The new analysis er- 
ror is equal to the background error reduced by a factor 
of (1 — KH). By finding the analysis, the filter has revealed 
the best possible starting point for the next background fore- 
cast. In fact, if the system is linear, the KF is the optimal 
algorithm for state-estimation. 

S2.2 Variational Filtering (3D-Var) 

Rather than minimize the analysis error variance, the 
analysis equations can also be derived by finding the analysis 
state x° that minimizes the quadratic scalar cost function 
2C(x) = (x - x'')^B-i(x - x*") + {y- UxfR-^y - Hx). 
The cost C(x) has its minimum at x = x" , where x" is given 
by Eqn. (fS^ . This is caUed the 3D variational (3D-Var) 
method since the minimization for NWP is with respect to a 
state vector embedded in a three-dimensional field (latitude, 
longitude, and height). 

Formally , both 3D-Var and the KF yield the same solu- 
tion (|Kalna4[200i l. However, in this case the control vari- 
able is the analysis, while in the KF the control variable 
is the weight matrix itself. In operational NWP, where the 
dimension of the state space iV is of 0(10^), the numerical 
implementations of 3D-Var and the nonlinear KF are dras- 
tically different. Because 3D-Var assumes the background 
error B is fixed in time, the Kalman gain K needs to be cal- 
culated only once. The calculation of K is the most compu- 
tationally prohibitive part of DA because it requires solving 



a linear system in A*" variables. A constant K thus makes the 
algorithm computationally simple; the most difficult part of 
implementing 3D-Var is finding the optimal B. 

However, a static B is not realistic. From a dynami- 
cal systems standpoint, uncertainty is closely related to sta- 
bility, which is clearly dependent on the system state. In 
the thermosyphon, the true background error is typically 
smaller when the system state is near the unstable convect- 
ing equilibria than when the state is near the more unstable 
conducting equilibrium. Because 3D-Var is computationally 
cheap, the National Centers for Environmental Prediction 
(NCEP) employ it to estimate ICs for the National Weather 
Service 14-day global forecasts. However, it cannot detect 
so-called "errors of the day" , state-dependent forecast errors 
which grow quickly but are not rep resented in the 3D-Var 
background error covariance matrix (|Kalnavl . [200^ :1 Li et al.l . 
l2009h . 

In our implementation, the 3D-Var background error 
covariance Bgo-Var was calculated i teratively, u s ing a tech- 
niques similar to that described in lYang et al.l (|2006l ') . We 
did this by calculating a time average of the outer product 
of analysis increments 

B3D-Var = {(x'' - x")(x'' - x")^), (S2-7) 

disregarding the initial 500 assimilation cycles and iterat- 
ing the process until convergence. During this, forecast er- 
rors were observed to decrease and stabilize. This BsD-Var 
was first computed for the 30 s assimilation window. It was 
stored and then used to bootstrap the iterative procedure 
for the 60 s assimilation window, which was stored and fed 
into the calculation for the 90 s assimilation window, etc. 



S2.3 Extended Kalman Filter (EKF) 

The EKF is essentially the KF applied to a nonlinear 
model. Given a nonlinear model M, the error covariances 
are updated by the linear tangent model M = dM /d'x.\^^^b 
which takes the place of M in Eqn. (|S2-4b[) . This model 
propagates small perturbations around the trajectory x*" for- 
ward in time. To operate on the matrix A with the linear 
tangent model, first take the Jacobian of F (the right hand 
side of the nonlinear differential equation x = -F(x) which 
describes the model A^) and evaluate it at the background 
point x*"; call this matrix J. Each column of A, which can 
be thought of as an error perturbation to the analysis state, 
is then integrated forward in time according to the linear 
ODE &i = Jan. 

Also note that if the observation operator H is nonlin- 
ear, it is replaced by a similar linear tangent model H in the 
matrix equations HS2-5fl and (|S2-6P . The transpose of these 
matrix functions are called adjoint models, which are used 
in sensitivity analysis of the state to perturbations. 

To propagate the background covariance without the 
explicit adjoint model, as Eqn. (|S2-4b|l would require, 
B was first decomposed w ith the Cholesky factorization 
(|Golub and van Loanl . Il996l ) into the product of a lower and 
upper diagonal matrix before its columns were integrated 
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forward with the hnear tangent model M. 

Tfc = Mfc_iLfe_i 



Q 



(S2-8) 
(S2-9) 
(S2-10) 



This guarantees symmetry for the new analysis error covari- 
ance A. 

Some modifications to the EKF algorithm are necessary 
to prevent filter divergence. A multiplicative inflation factor 



B 



(1 + A)B 



(S2-ff) 



was applied to the background covariance matrix after the 
model integration and before the a nalysis step. We al so per- 
formed additive inflation, following I Yang et al.l (|2006l V Ran- 
dom numbers uniformly distributed between and jj, were 
added to the diagonal elements of A after performing the 
analysis and before the next forecast step, i.e. 



A A + /i diag(!/) 



(S2-12) 



where u is an A^-dimensional vector whose entries are drawn 
from a uniform distribution between and 1. 



S2.4 Ensemble Kalman Filters (EnKFs) 

The EnKF is a method that replaces a single forecast 
state with an ensemble of states. The spread of the ensemble 
about its mean gives an approximation of the background 
error covariance and forecast uncertainty, while the ensem- 
ble average gives the be st guess o f the forecast. The EnKF 
was first introduced by lEvensenl (1 1994 ). For a com prehen- 
sive overview of ensemble filters, see lEvensenl (|2003h . It was 
shown that if the observation, which has random error with 
covariance R, is perturbed with P random errors (again with 
covariance R), to make an P-member ensemble of indepen- 
dent observations { yi|, then the b ackground error covari- 
ance can be written (lEvensenl. lioO^ ) 



B 



X'X'^ (S2-13) 



which is simply the unbiased average outer product of back- 



ground perturbations X 



The background 



forecast of ensemble member i is denoted , x*" is the back- 



ground forecast ensemble average, and = x^ 



x*" is the 



ith member's deviation from the mean. In this case, each 
ensemble member is updated according to the KF equations 
for their associated observation 

The advantages of the EnKF are many: there is no linear 
tangent model to compute, the number of ensemble mem- 
bers can be small (0(10'^) for NWP) relative to the dimen- 
sionality of the state space, and prior knowledge about the 
structure of the forecast errors is not necessary. Currently, 
4D-Var (like 3D-Var but also taking into account older ob- 
servations) and ensemble filters are the most promising can- 
didates being considered to replace 3D-Var in operational 
NWP. 

As with the EKF, ensemble filters tend to under- 
estimate the background error, resulting in an ensemble 
spread which is typically less too small. We again used 
multiplicative inflation of the backgroun d error, a com- 
mon method shown to be successful in lEvensenl (|2003D : 



IWhitaker and Hamilll ( 


20021): Annan and Hargreavesl 


(2004h: Yang et al.l (|200e 


): Kalnav et all (20071). This is 



accomplished by setting 



X'^ (1 + A)''''X 



(S2-14) 



before the analysis step. Additive inflation proved crucial to 
stabilizing both EnKFs tested. Without it, the filters some- 
times worked but only with A ^ 1; A is supposed to be 
a small parameter. As in the EKF, additive inflation is ap- 
plied immediately after the analysis step, but in this case 
the noise is added to the analysis ensemble states 



(S2-15) 



for all i = 1, . . . , P. The noise u is, again, an Ai'-dimensional 
random vector with entries drawn from the uniform distri- 
bution between and 1. 



S2.5 Ensemble Square Root Filter (EnSRF) 

The original EnKF adds noise to create linearly inde- 
pendent observ ations and is c l assifie d as a perturbed observa- 
tions method (|Kalnav et al.l . [2007h . This necessarily intro- 
duces additional sampling er r or int o the forecast. For this 
reason, IWhitaker and Hamilll (|2002h introduced the ensem- 
ble square root filter (EnSRF) as an improved EnKF. In the 
EnSRF, the ensemble mean is updated with the traditional 
Kalman gain (Eqn. ()S2-6|) ) 

5i^ = x^ + K(y--Hxfc) (S2-16) 

and deviations from the mean are updated by 



X" = (1 - KHW 



(S2-17) 



where 



K = BH' 



\/hBH^ + R 



X VHBH' +R + VR 



(S2-18) 

When the observation is a scalar, it can be shown that 



K 



l-f- 



R 



HBH" +R 



K 



(S2-19) 



If observation errors are uncorrelated (R is diagonal), then 
Eqn. (|S2-19ll can be used to proce ss observations one at 
a time fjWhitaker and Hamilll . |2002| ) . The updated analysis 
ensemble is then {x°}, where x" = x" -f x^". Square root 
filters have better numerical stability and speed than their 
standard KF counterparts. The Potter square root filter was 
employed for navigation in t he Lunar Module of the Apollo 
program (jSavelv et al.l.ll972h . 



S2.6 Ensemble Transform Kalman Filter (ETKF) 

The ETKF is another type of deterministic square root 
filter. In this variant, the analysis perturbations are assumed 
to be equal to the background perturbations postmultiplied 
by a transformation matrix T so that the analysis error co- 
variance satisfies Eqn. ()S2-5b|l . The analysis covariance is 
written 



1 



P- 1 



_^a^aT ^ x'AX'' 
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ensrf 360 s 




Figure S2-1. Contour plot of forecast error using the EnSRF and 
an assimilation window of 360 s, as the inflation parameters A 
and fj, were varied. The chosen parameters A = 0.15 and fi = 0.25 
lie in the region of lowest error. 

where A = [(P- 1)1+ (HX'')^R-i (HX'')]-^ . The analysis 
perturbations are X " = X^'T, where T = [(P- 1)A]^''^ See 
iKalnav et al] (|2007l 'l for further details. 

The local ensemble transform filter (LETKF) is a vari- 
ant that computes the analysis at a given gridpoint using 
only local observations. This allows for efficient paralleliza- 
tion. Localization removes spurious long-distance correla- 
tions from B and allows greater flexibility in the global anal- 
ysis by allowing diflerent linear combi nations of ensemble 
members at differ ent spatial locations (|Kalnav et al I l2007l: 
iHunt et al.l.l200'it ). 



S2.7 Tuning Parameters 

Table [S2^ lists the tuning parameters used for the DA 
experiments. The tuning was done manually. Sensitivity of 
model error to the tuning parameters was checked by creat- 
ing a course contour plot of background error for assimila- 
tion windows of 2, 4, 6, 8, and 10 minutes for each filter; an 
example is shown in Fig. IS2-1I 
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Table S2-1. Inflation tuning parameters used in DA experiments. 



APPENDIX S3: Flow reversal forecast tuning 

In Fig. IS3-1I we present skill score curves as the tuning 
parameters of the three flow reversal forecasts are varied. We 
used a sequence of plots of this type to inform our tuning of 
the various flow reversal tests. 

In Fig. IS3-2I we also present warning time histograms 
for the different tests. If early detection of flow reversals 
is desirable, then the warning times tell us how the tests 
compare. 
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Figure S3-1. Tuning tiic flow reversal forecasts metiiods by vary- 
ing (a) pcorr for tlie correlation test, (b) pBV for the BV test, and 
(c) Ajgad for lead test. The dashed line is at 95%. There are trade- 
offs among the various skill scores (POD, TS, RPS-avg, RPS- 
med) and costs (FAR). The chosen, final parameters Ai^ad = 7, 
Pbv = 0.6802, and pcorr = 1-4 (with Acorr = 18) appear in the 
middle of each plot. 
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Figure S3-2. Warning time distributions for the three flow re- 
versal tests, in units of 30 s analysis cycles. The warning time is 
the time interval between the test being triggered and observa- 
tion of the the actual flow reversal. The correlation test gives the 
earliest average warning time, followed by the BV and lead tests. 
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